# Annotation --------------------------------------------------------------
### This script replicates the analyses published in Toshkov, Mader and Rasmussen (2018), Journal of Public Policy
### 'Party Government and Policy Responsiveness. Evidence from Three Parliamentary Democracies'
### The analyses have been last run using R version 3.5.1 (2018-07-02) on 26 September 2018
# Set working directory ----------------------------------------------------------
setwd("E://PAPERS/ONGOING/Legislative Politics and Congruence (GOVLIS)/")
# Load libraries ----------------------------------------------------------
library(car) #data transformations
library(zoo) #date transformations
library(boot)#for the inverse logit function
library(survival) #survival analysis
library(effects) #illustration of interaction effects
library(gridExtra) #plotting
# Load data ---------------------------------------------------------------
# 1.Issue-level file with policy issue data
d<-read.table("./data/policy_issue_data.txt", sep='\t', head=T)
#recode as dates
d$po.date<-as.Date(d$po.date)
d$po.date.m<-as.yearmon(d$po.date.m)
d$policy.date<-as.Date(d$policy.date)
d$policy.date.m<-as.yearmon(d$policy.date.m)
d$policy.enact.date.m<-as.yearmon(d$policy.enact.date.m)
# 2. File containing the concrete issue-level party and government positions in Germany
new<-read.table("./data/party_positions_de_handcoded.txt", sep='\t', head=T)
# 3. File with Chapel Hill government positions data 
govpositions<-read.table("./data/govpositions2.txt", sep='\t', head=T) 
govpositions$edate<-as.Date(govpositions$edate)
govpositions$edate.m<-as.yearmon(govpositions$edate,"%b %Y")
govpositions$exp.end<-as.Date(govpositions$exp.end)
govpositions$exp.end.m<-as.yearmon(govpositions$exp.end,"%b %Y")
# Prepare a matrix with the month as a row and the country and government index as columns ----------------------------------
len<-length(seq(min(d$po.date.m),max(d$policy.date.m),1/12)) #get all months between the start of the earliest poll and the end of the last monitoring

#create an empty holder matrix
govs<-data.frame(matrix(ncol = 3, nrow =len ))
colnames(govs)<-c("month","country","govindex") #assign the column names  
govs$month<-seq(min(d$po.date.m),max(d$policy.date.m),1/12) #populate the month

#Populate the matrix for Germany
govs$country<-rep(c("DE"), each=len) #create a country column
# the rest assigns the government number
govs[1:(as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DE" & govpositions$govindex==2, "edate.m"],])))-1), "govindex"]<-1
govs[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DE" & govpositions$govindex==2, "edate.m"],]))):
       (as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DE" & govpositions$govindex==3, "edate.m"],])))-1), "govindex"]<-2
govs[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DE" & govpositions$govindex==3, "edate.m"],]))):
       (as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DE" & govpositions$govindex==4, "edate.m"],])))-1), "govindex"]<-3        
govs[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DE" & govpositions$govindex==4, "edate.m"],]))):
       (as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DE" & govpositions$govindex==5, "edate.m"],])))-1), "govindex"]<-4        
govs[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DE" & govpositions$govindex==5, "edate.m"],]))):
       (as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DE" & govpositions$govindex==6, "edate.m"],])))-1), "govindex"]<-5        
govs[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DE" & govpositions$govindex==6, "edate.m"],]))):
       len, "govindex"]<-6  

#now the same for UK
govs.uk<-data.frame(matrix(ncol = 3, nrow =len ))
colnames(govs.uk)<-c("month","country","govindex")   
govs.uk$month<-seq(min(d$po.date.m),max(d$policy.date.m),1/12)
govs.uk$country<-rep(c("UK"), each=len)
govs.uk[1:(as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="UK" & govpositions$govindex==2, "edate.m"],])))-1), "govindex"]<-1
govs.uk[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="UK" & govpositions$govindex==2, "edate.m"],]))):
          (as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="UK" & govpositions$govindex==3, "edate.m"],])))-1), "govindex"]<-2
govs.uk[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="UK" & govpositions$govindex==3, "edate.m"],]))):
          (as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="UK" & govpositions$govindex==4, "edate.m"],])))-1), "govindex"]<-3        
govs.uk[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="UK" & govpositions$govindex==4, "edate.m"],]))):
          len, "govindex"]<-4  

#and the same for Denmark
govs.dk<-data.frame(matrix(ncol = 3, nrow =len ))
colnames(govs.dk)<-c("month","country","govindex")   
govs.dk$month<-seq(min(d$po.date.m),max(d$policy.date.m),1/12)
govs.dk$country<-rep(c("DK"), each=len)
govs.dk[1:(as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==2, "edate.m"],])))-1), "govindex"]<-1
govs.dk[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==2, "edate.m"],]))):
          (as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==3, "edate.m"],])))-1), "govindex"]<-2
govs.dk[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==3, "edate.m"],]))):
          (as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==4, "edate.m"],])))-1), "govindex"]<-3        
govs.dk[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==4, "edate.m"],]))):
          (as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==5, "edate.m"],])))-1), "govindex"]<-4        
govs.dk[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==5, "edate.m"],]))):
          (as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==6, "edate.m"],])))-1), "govindex"]<-5        
govs.dk[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==6, "edate.m"],]))):
          (as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==7, "edate.m"],])))-1), "govindex"]<-6        
govs.dk[as.numeric(as.character(row.names(govs[govs$month==govpositions[govpositions$country=="DK" & govpositions$govindex==7, "edate.m"],]))):
          len, "govindex"]<-7  
#finally, combine the govpositions
govs.all<-rbind(govs, govs.dk, govs.uk)
# Now work this matrix to get new positions assigned per each case ----------------------------------
#create the additional varaibles
for ( i in 1:nrow(d)){ #for every case/ row in d
  temp.m<-length(seq(d$po.date.m[i],d$policy.date.m[i],1/12)) #get the number of months between begining and end
  temp<-as.data.frame(matrix(NA,temp.m, 33)) #create a matrix with that amount of rows and the required number of columns
  temp[,1]<-seq(d$po.date.m[i],d$policy.date.m[i],1/12) #put the month first
  temp[,2]<-d[i, "country"] #then the country
  for (j in 1:nrow(temp)){ #then for each month (row in the temporary matrix)
    temp[j,3]<-govs.all[govs.all$country==temp[j,2] & govs.all$month==temp[j,1],"govindex"] #assign the government index number
  }
  for (k in 1:nrow(temp)){ #now put the relevant variables getting them from
    temp[k,4]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"minlrgen"]
    temp[k,5]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"maxlrgen"]
    temp[k,6]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"govposlrgen"]
    temp[k,7]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"minlrecon"]
    temp[k,8]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"maxlrecon"]
    temp[k,9]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"govposlrecon"]
    temp[k,10]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"mingaltan"]
    temp[k,11]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"maxgaltan"]
    temp[k,12]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"govposgaltan"]
    
    temp[k,13]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"minlrgenDK"]
    temp[k,14]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"maxlrgenDK"]
    temp[k,15]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"govposlrgenDK"]
    temp[k,16]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"minlreconDK"]
    temp[k,17]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"maxlreconDK"]
    temp[k,18]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"govposlreconDK"]
    temp[k,19]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"mingaltanDK"]
    temp[k,20]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"maxgaltanDK"]
    temp[k,21]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"govposgaltanDK"]
    
    temp[k,22]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"n.parties"]
    temp[k,23]<-govpositions[govpositions$country==temp[k,2] & govpositions$govindex==temp[k,3],"n.partiesDK"]
    
    if (temp[k,2]=='DE' & temp[k,3]==1)
    {temp[k,30]<-new[new$id==d$IDde[i], "mingov1"]
    temp[k,31]<-new[new$id==d$IDde[i], "maxgov1"]  
    temp[k,32]<-new[new$id==d$IDde[i], "meangov1"]
    temp[k,33]<-new[new$id==d$IDde[i], "range1"]}
    else
      if (temp[k,2]=='DE' & temp[k,3]==2)
      {temp[k,30]<-new[new$id==d$IDde[i], "mingov2"]
      temp[k,31]<-new[new$id==d$IDde[i], "maxgov2"]  
      temp[k,32]<-new[new$id==d$IDde[i], "meangov2"]
      temp[k,33]<-new[new$id==d$IDde[i], "range2"]}
    else
      if (temp[k,2]=='DE' & temp[k,3]==3)
      {temp[k,30]<-new[new$id==d$IDde[i], "mingov3"]
      temp[k,31]<-new[new$id==d$IDde[i], "maxgov3"]  
      temp[k,32]<-new[new$id==d$IDde[i], "meangov3"]
      temp[k,33]<-new[new$id==d$IDde[i], "range3"]}
    else
      if (temp[k,2]=='DE' & temp[k,3]==4)
      {temp[k,30]<-new[new$id==d$IDde[i], "mingov4"]
      temp[k,31]<-new[new$id==d$IDde[i], "maxgov4"]  
      temp[k,32]<-new[new$id==d$IDde[i], "meangov4"]
      temp[k,33]<-new[new$id==d$IDde[i], "range4"]}
    else
      if (temp[k,2]=='DE' & temp[k,3]==5)
      {temp[k,30]<-new[new$id==d$IDde[i], "mingov5"]
      temp[k,31]<-new[new$id==d$IDde[i], "maxgov5"]  
      temp[k,32]<-new[new$id==d$IDde[i], "meangov5"]
      temp[k,33]<-new[new$id==d$IDde[i], "range5"]}
    else
      if (temp[k,2]=='DE' & temp[k,3]==6)
      {temp[k,30]<-new[new$id==d$IDde[i], "mingov6"]
      temp[k,31]<-new[new$id==d$IDde[i], "maxgov6"]  
      temp[k,32]<-new[new$id==d$IDde[i], "meangov6"]
      temp[k,33]<-new[new$id==d$IDde[i], "range6"]}
    else
    {temp[k,30]<-NA
    temp[k,31]<-NA
    temp[k,32]<-NA
    temp[k,33]<-NA}
  }
  
  temp[,24]<-abs(temp[,5]-temp[,4]) #get the ranges
  temp[,25]<-abs(temp[,8]-temp[,7])
  temp[,26]<-abs(temp[,11]-temp[,10])
  
  temp[,27]<-abs(temp[,14]-temp[,13]) #get the ranges
  temp[,28]<-abs(temp[,17]-temp[,16])
  temp[,29]<-abs(temp[,20]-temp[,19])
  
  #now assign the mean of the rows of the temporary matrix as entries in the general dataset
  d$minlrgen[i]<-mean(temp[,4])
  d$maxlrgen[i]<-mean(temp[,5])
  d$govposlrgen[i]<-mean(temp[,6])
  d$rangelrgen[i]<-mean(temp[,24])
  
  d$minlrecon[i]<-mean(temp[,7])
  d$maxlrecon[i]<-mean(temp[,8])
  d$govposlrecon[i]<-mean(temp[,9])
  d$rangelrecon[i]<-mean(temp[,25])
  
  d$mingaltan[i]<-mean(temp[,10])
  d$maxgaltan[i]<-mean(temp[,11])
  d$govposgaltan[i]<-mean(temp[,12])
  d$rangegaltan[i]<-mean(temp[,26])
  
  d$averageparties[i]<-mean(temp[,22])
  
  d$minlrgenDK[i]<-mean(temp[,13])
  d$maxlrgenDK[i]<-mean(temp[,14])
  d$govposlrgenDK[i]<-mean(temp[,15])
  d$rangelrgenDK[i]<-mean(temp[,27])
  
  d$minlreconDK[i]<-mean(temp[,16])
  d$maxlreconDK[i]<-mean(temp[,17])
  d$govposlreconDK[i]<-mean(temp[,18])
  d$rangelreconDK[i]<-mean(temp[,28])
  
  d$mingaltanDK[i]<-mean(temp[,19])
  d$maxgaltanDK[i]<-mean(temp[,20])
  d$govposgaltanDK[i]<-mean(temp[,21])
  d$rangegaltanDK[i]<-mean(temp[,29])
  
  d$averagepartiesDK[i]<-mean(temp[,23])
  
  d$real_min[i]<-mean(temp[,30], na.tm=T)
  d$real_max[i]<-mean(temp[,31], na.tm=T)
  d$real_mean[i]<-mean(temp[,32], na.tm=T)
  d$real_range[i]<-mean(temp[,33], na.tm=T)
}

#now assign the numbers from the relevant policy scale
d$govpos.main<-ifelse(d$policy.scale==1, d$govposlrgen, ifelse(d$policy.scale==2, d$govposlrecon, d$govposgaltan)) 
d$min.main<-ifelse(d$policy.scale==1, d$minlrgen, ifelse(d$policy.scale==2, d$minlrecon, d$mingaltan)) 
d$max.main<-ifelse(d$policy.scale==1, d$maxlrgen, ifelse(d$policy.scale==2, d$maxlrecon, d$maxgaltan)) 
d$range.main<-ifelse(d$policy.scale==1, d$rangelrgen, ifelse(d$policy.scale==2, d$rangelrecon, d$rangegaltan)) 
d$govpos.main.c<-d$govpos.main-5 #center
d$min.main.c<-d$min.main-5 #center
d$max.main.c<-d$max.main-5 #center

#create the goverment support for specific policies variables
#for right-going policies, this is the average government positions AND the minimum of the partners
#for left-going policies, this is the negative of the average government position AND the negative of the maximum of the partners

d$gov.support<-ifelse(d$policy.direction.right==1, d$govpos.main.c, -d$govpos.main.c)
d$gov.vetosupport<-ifelse(d$policy.direction.right==1, d$min.main.c, -d$max.main.c)

d$govpos.mainDK<-ifelse(d$policy.scale==1, d$govposlrgenDK, ifelse(d$policy.scale==2, d$govposlreconDK, d$govposgaltanDK)) 
d$min.mainDK<-ifelse(d$policy.scale==1, d$minlrgenDK, ifelse(d$policy.scale==2, d$minlreconDK, d$mingaltanDK)) 
d$max.mainDK<-ifelse(d$policy.scale==1, d$maxlrgenDK, ifelse(d$policy.scale==2, d$maxlreconDK, d$maxgaltanDK)) 
d$range.mainDK<-ifelse(d$policy.scale==1, d$rangelrgenDK, ifelse(d$policy.scale==2, d$rangelreconDK, d$rangegaltanDK)) 
d$govpos.mainDK.c<-d$govpos.mainDK-5 #center
d$min.mainDK.c<-d$min.mainDK-5 #center
d$max.mainDK.c<-d$max.mainDK-5 #center

d$gov.supportDK<-ifelse(d$policy.direction.right==1, d$govpos.mainDK.c, -d$govpos.mainDK.c)
d$gov.vetosupportDK<-ifelse(d$policy.direction.right==1, d$min.mainDK.c, -d$max.mainDK.c)


# Create a dataset t with month as the row ----------------------------------
d$country<-as.character(d$country) #to preserve the levels
nc<-ncol(d)
t<-as.data.frame(matrix(NA,0, nc+3))
colnames(t)<-c(colnames(d),"case.n","month.c","month")
for (i in 1:nrow(d)){
  n<-max((d$policy.date.m[i]-d$po.date.m[i])*12,1)
  last<-nrow(t)+1
  t[last:(last+n-1),1:nc]<-d[i,]
  t[last:(last+n-1),"case.n"]<-paste("case", i, sep=".")
  t[last:(last+n-1),"month.c"]<-1:n
  t[last:(last+n-1),"month"]<-as.Date(d[i,"po.date.m"]+(1/12)*(1:n))
}


#reconstruct the time of public opinion as a date and as a month
t$po.date<-as.Date(t$po.date, origin="1970-01-01")
t$po.date.m<-as.yearmon(t$po.date,"%b %Y")
t$po.expiry.m<-t$po.date.m+4
# Add the number of government --------------------------------------------
for (i in 1:nrow(t)){
  temp<-subset(govpositions, country==t$country[i])
  l<-nrow(temp)
  if (t$country[i]=="DE") {
    if(as.Date(t$month[i])<govpositions[govpositions$country=="DE" & govpositions$govindex==2,"edate"])
      t$government[i]<-1
    else
      if(as.Date(t$month[i])<govpositions[govpositions$country=="DE" & govpositions$govindex==3,"edate"])
        t$government[i]<-2
      else
        if(as.Date(t$month[i])<govpositions[govpositions$country=="DE" & govpositions$govindex==4,"edate"])
          t$government[i]<-3
        else
          if(as.Date(t$month[i])<govpositions[govpositions$country=="DE" & govpositions$govindex==5,"edate"])
            t$government[i]<-4
          else
            if(as.Date(t$month[i])<govpositions[govpositions$country=="DE" & govpositions$govindex==6,"edate"])
              t$government[i]<-5
            else
              t$government[i]<-6}
  else
    if (t$country[i]=="DK") {
      if(as.Date(t$month[i])<govpositions[govpositions$country=="DK" & govpositions$govindex==2,"edate"])
        t$government[i]<-1
      else
        if(as.Date(t$month[i])<govpositions[govpositions$country=="DK" & govpositions$govindex==3,"edate"])
          t$government[i]<-2
        else
          if(as.Date(t$month[i])<govpositions[govpositions$country=="DK" & govpositions$govindex==4,"edate"])
            t$government[i]<-3
          else
            if(as.Date(t$month[i])<govpositions[govpositions$country=="DK" & govpositions$govindex==5,"edate"])
              t$government[i]<-4
            else
              if(as.Date(t$month[i])<govpositions[govpositions$country=="DK" & govpositions$govindex==6,"edate"])
                t$government[i]<-5
              else
                if(as.Date(t$month[i])<govpositions[govpositions$country=="DK" & govpositions$govindex==7,"edate"])
                  t$government[i]<-6
                else
                  t$government[i]<-7
    }
  else
    if (t$country[i]=="UK") {
      if(as.Date(t$month[i])<govpositions[govpositions$country=="UK" & govpositions$govindex==2,"edate"])
        t$government[i]<-1
      else
        if(as.Date(t$month[i])<govpositions[govpositions$country=="UK" & govpositions$govindex==3,"edate"])
          t$government[i]<-2
        else
          if(as.Date(t$month[i])<govpositions[govpositions$country=="UK" & govpositions$govindex==4,"edate"])
            t$government[i]<-3
          else
            t$government[i]<-4
    }
  
}

t$exp.end<-t$po.date #do this first  to preserve the date format of the variables created below 
t$exp.end.m<-t$po.expiry.m #do this first  to preserve the date format of the variables created below 
t$edate<-t$po.date #do this first  to preserve the date format of the variables created below 
t$edate.m<-t$po.expiry.m #do this first  to preserve the date format of the variables created below 

#Now assign the actual expected termination dates of governments
for (i in 1:nrow(t)){
  t$exp.end[i]<-govpositions[govpositions$country==t[i,"country"] & govpositions$govindex==t[i,"government"],"exp.end"]
  t$exp.end.m[i]<-govpositions[govpositions$country==t[i,"country"] & govpositions$govindex==t[i,"government"],"exp.end.m"]
  t$edate[i]<-govpositions[govpositions$country==t[i,"country"] & govpositions$govindex==t[i,"government"],"edate"]
  t$edate.m[i]<-govpositions[govpositions$country==t[i,"country"] & govpositions$govindex==t[i,"government"],"edate.m"]
}

t$month.f<-as.Date(as.Date(t$month),"%b %Y")  
# Get the position variables --------------------------------------------
for (i in 1:nrow(t)){
  t$min.lrgen[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"minlrgen"]
  t$max.lrgen[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"maxlrgen"]
  t$mean.lrgen[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"govposlrgen"]
  
  t$min.lrecon[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"minlrecon"]
  t$max.lrecon[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"maxlrecon"]
  t$mean.lrecon[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"govposlrecon"]
  
  t$min.galtan[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"mingaltan"]
  t$max.galtan[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"maxgaltan"]
  t$mean.galtan[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"govposgaltan"]
  
  t$n.parties[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"n.parties"]
  
  t$min.lrgen.dk[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"minlrgenDK"]
  t$max.lrgen.dk[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"maxlrgenDK"]
  t$mean.lrgen.dk[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"govposlrgenDK"]
  
  t$min.lrecon.dk[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"minlreconDK"]
  t$max.lrecon.dk[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"maxlreconDK"]
  t$mean.lrecon.dk[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"govposlreconDK"]
  
  t$min.galtan.dk[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"mingaltanDK"]
  t$max.galtan.dk[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"maxgaltanDK"]
  t$mean.galtan.dk[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"govposgaltanDK"]
  
  t$n.partiesDK[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"n.partiesDK"]
  
  #########
  if (t$country[i]=='DE' & t$government[i]==1)
  {t$real_min_t[i]<-new[new$id==t$IDde[i], "mingov1"]
  t$real_max_t[i]<-new[new$id==t$IDde[i], "maxgov1"]  
  t$real_mean_t[i]<-new[new$id==t$IDde[i], "meangov1"]
  t$real_range_t[i]<-new[new$id==t$IDde[i], "range1"]}
  else
    if (t$country[i]=='DE' & t$government[i]==2)
    {t$real_min_t[i]<-new[new$id==t$IDde[i], "mingov2"]
    t$real_max_t[i]<-new[new$id==t$IDde[i], "maxgov2"]  
    t$real_mean_t[i]<-new[new$id==t$IDde[i], "meangov2"]
    t$real_range_t[i]<-new[new$id==t$IDde[i], "range2"]}
  else
    if (t$country[i]=='DE' & t$government[i]==3)
    {t$real_min_t[i]<-new[new$id==t$IDde[i], "mingov3"]
    t$real_max_t[i]<-new[new$id==t$IDde[i], "maxgov3"]  
    t$real_mean_t[i]<-new[new$id==t$IDde[i], "meangov3"]
    t$real_range_t[i]<-new[new$id==t$IDde[i], "range3"]}
  else
    if (t$country[i]=='DE' & t$government[i]==4)
    {t$real_min_t[i]<-new[new$id==t$IDde[i], "mingov4"]
    t$real_max_t[i]<-new[new$id==t$IDde[i], "maxgov4"]  
    t$real_mean_t[i]<-new[new$id==t$IDde[i], "meangov4"]
    t$real_range_t[i]<-new[new$id==t$IDde[i], "range4"]}
  else
    if (t$country[i]=='DE' & t$government[i]==5)
    {t$real_min_t[i]<-new[new$id==t$IDde[i], "mingov5"]
    t$real_max_t[i]<-new[new$id==t$IDde[i], "maxgov5"]  
    t$real_mean_t[i]<-new[new$id==t$IDde[i], "meangov5"]
    t$real_range_t[i]<-new[new$id==t$IDde[i], "range5"]}
  else
    if (t$country[i]=='DE' & t$government[i]==6)
    {t$real_min_t[i]<-new[new$id==t$IDde[i], "mingov6"]
    t$real_max_t[i]<-new[new$id==t$IDde[i], "maxgov6"]  
    t$real_mean_t[i]<-new[new$id==t$IDde[i], "meangov6"]
    t$real_range_t[i]<-new[new$id==t$IDde[i], "range6"]}
  else
  {t$real_min_t[i]<-NA
  t$real_max_t[i]<-NA
  t$real_mean_t[i]<-NA
  t$real_range_t[i]<-NA}
}

#assign according to scale
for (i in 1:nrow(t)){
  if (t$policy.scale[i]==1){
    t$min.main[i]<-t$min.lrgen[i]
    t$max.main[i]<-t$max.lrgen[i]
    t$mean.main[i]<-t$mean.lrgen[i]
    t$coal.conflict.main[i]<-abs(t$max.lrgen[i]-t$min.lrgen[i])
    
    t$min.main.dk[i]<-t$min.lrgen.dk[i]
    t$max.main.dk[i]<-t$max.lrgen.dk[i]
    t$mean.main.dk[i]<-t$mean.lrgen.dk[i]
    t$coal.conflict.main.dk[i]<-abs(t$max.lrgen.dk[i]-t$min.lrgen.dk[i])
  }
  else
    if (t$policy.scale[i]==2){
      t$min.main[i]<-t$min.lrecon[i]
      t$max.main[i]<-t$max.lrecon[i]
      t$mean.main[i]<-t$mean.lrecon[i]
      t$coal.conflict.main[i]<-abs(t$max.lrecon[i]-t$min.lrecon[i])
      
      t$min.main.dk[i]<-t$min.lrecon.dk[i]
      t$max.main.dk[i]<-t$max.lrecon.dk[i]
      t$mean.main.dk[i]<-t$mean.lrecon.dk[i]
      t$coal.conflict.main.dk[i]<-abs(t$max.lrecon.dk[i]-t$min.lrecon.dk[i])
    }
  else
    if (t$policy.scale[i]==3){
      t$min.main[i]<-t$min.galtan[i]
      t$max.main[i]<-t$max.galtan[i]
      t$mean.main[i]<-t$mean.galtan[i]
      t$coal.conflict.main[i]<-abs(t$max.galtan[i]-t$min.galtan[i])
      
      t$min.main.dk[i]<-t$min.galtan.dk[i]
      t$max.main.dk[i]<-t$max.galtan.dk[i]
      t$mean.main.dk[i]<-t$mean.galtan.dk[i]
      t$coal.conflict.main.dk[i]<-abs(t$max.galtan.dk[i]-t$min.galtan.dk[i])
    }
  else
    next
}

t$event<-0
u<-unique(t$case.n)

for (i in 1:length(u)){
  if (t[t$case.n==u[i],][nrow(t[t$case.n==u[i],]),"month.c"]>=48)
    next
  else
    t[t$case.n==u[i],][nrow(t[t$case.n==u[i],]),"event"]<-1
}
#works with the exception of case 119 for some reaoson

#some additional variables
t$mean.main.c<-t$mean.main-5 #center
t$mean.main.dk.c<-t$mean.main.dk-5 #center
t$min.main.c<-t$min.main-5 #center
t$min.main.dk.c<-t$min.main.dk-5 #center
t$max.main.c<-t$max.main-5 #center
t$max.main.dk.c<-t$max.main.dk-5 #center

t$gov.support.t<-ifelse(t$policy.direction.right==1, t$mean.main.c, -t$mean.main.c)
t$gov.support.dk.t<-ifelse(t$policy.direction.right==1, t$mean.main.dk.c, -t$mean.main.dk.c)

t$gov.vetosupport.t<-ifelse(t$policy.direction.right==1, t$min.main.c, -t$max.main.c)
t$gov.vetosupport.dk.t<-ifelse(t$policy.direction.right==1, t$min.main.dk.c, -t$max.main.dk.c)

t$final.exp<-t$exp.end.m #just to keep the date format
t$start.per<-t$exp.end.m #just to keep the date format
t$final.exp<-pmin(t$exp.end.m, t$po.expiry.m) #get the minimum end date of the government from the end of term or end of observations
t$start.per<-pmax(t$edate.m, t$po.date.m) #get the maximum start date of the government from the government start or the public opinion survey

t$remaining<-(t$final.exp-t$start.per)*12 #calculate the remaining time of the government for the case at the start of the government
t$index<-paste(t$case.n, t$country, t$government, sep=".")


# Prepare a dataset with issue/government as the unit of analysis, based on 't' ---------------------------------------------------------------------
###
ti<-unique(t$index)
g<-as.data.frame(matrix(NA, nrow=length(ti),ncol=30))
for (i in 1:length(ti)){
  t2<-subset(t, t$index==ti[i])
  g[i,1]<-t2$index[1]
  g[i,2]<-t2$country[1]
  g[i,3]<-t2$case.n[1]
  g[i,4]<-max(t2$event)
  
  g[i,5]<-t2$po.support[1]
  g[i,6]<-t2$net.support[1]
  g[i,7]<-t2$po.support.alt[1]
  
  g[i,8]<-t2$mean.main[1]
  g[i,9]<-t2$mean.main.dk[1]
  g[i,10]<-t2$min.main[1]
  g[i,11]<-t2$min.main.dk[1]
  g[i,12]<-t2$max.main[1]
  g[i,13]<-t2$max.main.dk[1]
  g[i,14]<-t2$coal.conflict.main[1]
  g[i,15]<-t2$coal.conflict.main.dk[1]
  g[i,16]<-t2$n.parties[1]
  g[i,17]<-t2$n.partiesDK[1]
  
  g[i,18]<-t2$gov.support.t[1]
  g[i,19]<-t2$gov.support.dk.t[1]
  g[i,20]<-t2$gov.vetosupport.t[1]
  g[i,21]<-t2$gov.vetosupport.dk.t[1]
  
  g[i,22]<-t2$salience.init[1]
  g[i,23]<-t2$proposal[1]
  g[i,24]<-t2$policy.direction.right[1]
  g[i,25]<-t2$policy.scale[1]
  g[i,26]<-t2$remaining[1]
  
  g[i,27]<-t2$real_min_t[1]
  g[i,28]<-t2$real_max_t[1]
  g[i,29]<-t2$real_mean_t[1]
  g[i,30]<-t2$real_range_t[1]
  
}

g<-as.data.frame(g, stringsAsFactors = FALSE)
colnames(g)<-c("index","country", "case.n", "event", 
               "po.support", "net.support", "po.support.alt",
               
               "govpos.main","govpos.main.dk","min.main","min.main.dk","max.main","max.main.dk",
               "coal.conflict.main", "coal.conflict.main.dk","n.parties","n.parties.dk", 
               
               "gov.support","gov.support.dk", "gov.vetosupport","gov.vetosupport.dk", 
               
               "salience.init", "proposal", "policy.direction.right", "policy.scale", "remaining",
               "real_min_g","real_max_g","real_mean_g","real_range_g")
# Country-level and scale subsets of the data --------------------------------------------
d.uk<-subset(d, country=="UK")
d.de<-subset(d, country=="DE")
d.dk<-subset(d, country=="DK")
d.gen<-subset(d, policy.scale==1)
d.econ<-subset(d, policy.scale==2)
d.galtan<-subset(d, policy.scale==3)

g.uk<-subset(g, country=="UK")
g.de<-subset(g, country=="DE")
g.dk<-subset(g, country=="DK")
g.gen<-subset(g, policy.scale==1)
g.econ<-subset(g, policy.scale==2)
g.galtan<-subset(g, policy.scale==3)

# Descriptives (Supplementary Material A1) ------------------------------------------------------------
d.list=c("po.support", "net.support", "salience.init",
         "policy.change.duration.m")
tableA1a<-matrix(NA, ncol=6, nrow=length(d.list))
for (i in 1:length(d.list)){
  tableA1a[i,]<-c(d.list[i], round(c(min(d[d.list[i]],na.rm=T),mean(unlist(d[d.list[i]]),na.rm=T),
                                     median(unlist(d[d.list[i]]),na.rm=T),max(unlist(d[d.list[i]]),na.rm=T),
                                     sd(unlist(d[d.list[i]]),na.rm=T)), digits=2))
}
tableA1a
g.list=c("govpos.main", "min.main", "max.main")
tableA1b<-matrix(NA, ncol=6, nrow=length(g.list))
for (i in 1:length(g.list)){
  tableA1b[i,]<-c(g.list[i], round(c(min(g[g.list[i]],na.rm=T)-5,mean(unlist(g[g.list[i]]),na.rm=T)-5,
                                     median(unlist(g[g.list[i]]),na.rm=T)-5,max(unlist(g[g.list[i]]),na.rm=T)-5,
                                     sd(unlist(g[g.list[i]]),na.rm=T)), digits=2))
}
tableA1b
g.list=c("coal.conflict.main")
tableA1c<-matrix(NA, ncol=6, nrow=length(g.list))
for (i in 1:length(g.list)){
  tableA1c[i,]<-c(g.list[i], round(c(min(g[g.list[i]],na.rm=T),mean(unlist(g[g.list[i]]),na.rm=T),
                                     median(unlist(g[g.list[i]]),na.rm=T),max(unlist(g[g.list[i]]),na.rm=T),
                                     sd(unlist(g[g.list[i]]),na.rm=T)), digits=2))
}
tableA1c

round(prop.table(table(d$proposal)), digits=2)
round(prop.table(table(d$policy.change)), digits=2)
round(prop.table(table(d$policy.direction.right)), digits=2)
round(prop.table(table(d$policy.scale)), digits=2)

# Country plots I: Public opinion and the likelihood of policy change [Figure A3 top] --------------------------------------------
summary(glm(policy.change~po.support, data=d.uk, family='binomial'))
summary(glm(policy.change~po.support, data=d.de, family='binomial'))
summary(glm(policy.change~po.support, data=d.dk, family='binomial'))

#pdf("./plots/pdf/POandpolicy.pdf",width=9,height=9)
#png("./plots/png/POandpolicy.png",width=1300, height=1300, res=160)

#with public support
par(mfrow=c(2,2))
par(mar=c(5,5,2,1))
x<-seq(0,1,0.01)
y.uk<-inv.logit(-2.1517+x*1.8136)
y.de<-inv.logit(-0.5334+x*0.6854)
y.dk<-inv.logit(-0.6651+x*1.1193)
plot(x=x, y=y.uk, type="l", xlim=c(0,1), ylim=c(0,1), lty=3, lwd=2.5, col="blue", 
     xlab="Public support for policy change", ylab="Probability of policy change enacted within 4 years")
abline(h=0.5, lty=2, col="grey")
abline(v=0.5, lty=2, col="grey")
lines(x=x, y=y.de, lty=1, lwd=2.5, col="black")
lines(x=x, y=y.dk, lty=4, lwd=2.5, col="red")
text(x=0.10, y=0.45, "Germany", col="black")
text(x=0.10, y=0.2, "UK", col="blue")
text(x=0.8, y=0.63, "Denmark", col="red")
rug(d.de$po.support, col="black")
rug(d.uk$po.support, col="blue", side=1)
rug(d.dk$po.support, col="red", side=1)
x<-seq(-1,1,0.01)
y.uk<-inv.logit(-1.3774+x*0.9566)
y.de<-inv.logit(-0.2113+x*0.3446)
y.dk<-inv.logit(-0.1868+x*0.7359)

plot(x=x, y=y.uk, type="l", xlim=c(-1,1), ylim=c(0,1), lty=3, lwd=2.5, col="blue", 
     xlab="Net public support for policy change", ylab="Probability of policy change enacted within 4 years")
abline(h=0.5, lty=2, col="grey")
abline(v=0, lty=2, col="grey")
lines(x=x, y=y.de, lty=1, lwd=2.5, col="black")
lines(x=x, y=y.dk, lty=4, lwd=2.5, col="red")
text(x=-0.8, y=0.45, "Germany", col="black")
text(x=0.8, y=0.2, "UK", col="blue")
text(x=0.8, y=0.7, "Denmark", col="red")
rug(d.de$net.support, col="black")
rug(d.uk$net.support, col="blue", side=1)
rug(d.dk$net.support, col="red", side=1)
table(d$policy.change, d$po.support.dummy)
table(d$policy.change, d$net.support.dummy)
# Country plots II: Public opinion and congruence [Figure A3 bottom] --------------------------------------------
x<-0:4
y.de<-c(prop.table(table(d.de$congruent0))[2],
        prop.table(table(d.de$congruent1))[2],
        prop.table(table(d.de$congruent2))[2],
        prop.table(table(d.de$congruent3))[2],
        prop.table(table(d.de$congruent4))[2])
y.dk<-c(prop.table(table(d.dk$congruent0))[2],
        prop.table(table(d.dk$congruent1))[2],
        prop.table(table(d.dk$congruent2))[2],
        prop.table(table(d.dk$congruent3))[2],
        prop.table(table(d.dk$congruent4))[2])
y.uk<-c(prop.table(table(d.uk$congruent0))[2],
        prop.table(table(d.uk$congruent1))[2],
        prop.table(table(d.uk$congruent2))[2],
        prop.table(table(d.uk$congruent3))[2],
        prop.table(table(d.uk$congruent4))[2])

#par(mfrow=c(1,2))
plot(x=x, y=y.uk, type="p", pch=16, cex=1.5, xlim=c(0,4), ylim=c(0,1), lty=3, lwd=2.5, col="blue", xaxt="n",
     xlab="Congruence between majority opinion and policy over time", ylab="Share of policies congruent with majority public opinion")
abline(h=seq(0,1,0.25), lty=2, col="grey")
points(x=x, y=y.de, pch=15, cex=1.5, col="black")
points(x=x, y=y.dk, pch=17, cex=1.5, col="red")
lines(x=x, y=y.de, lty=1, col="black")
lines(x=x, y=y.uk, lty=3, col="blue")
lines(x=x, y=y.dk, lty=4, col="red")
text(x=3.5, y=0.55, "Germany", col="black")
text(x=3.5, y=0.45, "UK", col="blue")
text(x=3.5, y=0.65, "Denmark", col="red")
axis(1, at=0:4, labels=c("At start","+1 year","+2 years","+3 years","+4 years"))

x<-0:4
y.de<-c(0, prop.table(table(d.de$policy.change.1))[2],
        prop.table(table(d.de$policy.change.2))[2],
        prop.table(table(d.de$policy.change.3))[2],
        prop.table(table(d.de$policy.change))[2])
y.uk<-c(0, prop.table(table(d.uk$policy.change.1))[2],
        prop.table(table(d.uk$policy.change.2))[2],
        prop.table(table(d.uk$policy.change.3))[2],
        prop.table(table(d.uk$policy.change))[2])

y.dk<-c(0, prop.table(table(d.dk$policy.change.1))[2],
        prop.table(table(d.dk$policy.change.2))[2],
        prop.table(table(d.dk$policy.change.3))[2],
        prop.table(table(d.dk$policy.change))[2])

plot(x=x, y=y.uk, type="p", pch=16, cex=1.5, xlim=c(0,4), ylim=c(0,1), lwd=2.5, col="blue", xaxt="n",
     xlab="Policy change over time", ylab="Share of public policies changed over time")
abline(h=seq(0,1,0.25), lty=2, col="grey")
points(x=x, y=y.de, pch=15, cex=1.5, col="black")
points(x=x, y=y.dk, pch=17, cex=1.5, col="red")
lines(x=x, y=y.de, lty=1, col="black")
lines(x=x, y=y.uk, lty=3, col="blue")
lines(x=x, y=y.dk, lty=4, col="red")
text(x=3.5, y=0.35, "Germany", col="black")
text(x=3.5, y=0.15, "UK", col="blue")
text(x=3.5, y=0.55, "Denmark", col="red")
axis(1, at=0:4, labels=c("At start","+1 year","+2 years","+3 years","+4 years"))
#dev.off()

# Plot of public opinion effects with lowess [Figure A4] --------------------------------------------
#pdf("./plots/pdf/POloess.pdf",width=9,height=7)
#png("./plots/png/POloess.png",width=1300, height=800, res=120)

par(mfrow=c(1,2))
par(mar=c(5,5,2,1))
#A. with public support
plx<-predict(loess(d$policy.change~d$po.support), se=T)
plx.up<-plx$fit + qt(0.975,plx$df)*plx$se
plx.down<-plx$fit - qt(0.975,plx$df)*plx$se
temp<-cbind(d$po.support,plx$fit, plx.up, plx.down)
temp<-temp[order(temp[,1]),]
plot(temp[,1], temp[,2], type="l", lty=1, lwd=2.5, xlab="Public support", ylab="Predicted policy change (loess)", ylim=c(0,1), xlim=c(0,1), col="black")
lines(temp[,1], temp[,3], type="l", lty=2, lwd=1.5)
lines(temp[,1], temp[,4], type="l", lty=2, lwd=1.5)
points(d.uk$po.support,predict(loess(d.uk$policy.change~d.uk$po.support)), col="lightblue", cex=0.8)
points(d.dk$po.support,predict(loess(d.dk$policy.change~d.dk$po.support)), col="coral", cex=0.8)
points(d.de$po.support,predict(loess(d.de$policy.change~d.de$po.support)), col="grey", cex=0.8)
lines(temp[,1], temp[,2], type="l", lty=1, lwd=3.5)
text(x=0.6, y=0.1, "UK", col='lightblue')
text(x=0.3, y=0.6, "Germany", col='grey')
text(x=0.7, y=0.7, "Denmark", col='coral')
text(x=0.5, y=0.30, "All cases", col='black')

#B. and with net suppport
plx2<-predict(loess(d$policy.change~d$net.support), se=T)
plx2.up<-plx2$fit + qt(0.975,plx2$df)*plx2$se
plx2.down<-plx2$fit - qt(0.975,plx2$df)*plx2$se
temp2<-cbind(d$net.support,plx2$fit, plx2.up, plx2.down)
temp2<-temp2[order(temp2[,1]),]
plot(temp2[,1], temp2[,2], type="l", lwd=2.5, xlab="Net public support", ylab="Predicted policy change (loess)", ylim=c(0,1), xlim=c(-1,1), col="black")
lines(temp2[,1], temp2[,3], type="l", lty=2, lwd=1.5)
lines(temp2[,1], temp2[,4], type="l", lty=2, lwd=1.5)
points(d.uk$net.support,predict(loess(d.uk$policy.change~d.uk$net.support)), col="lightblue", cex=0.8)
points(d.dk$net.support,predict(loess(d.dk$policy.change~d.dk$net.support)), col="coral", cex=0.8)
points(d.de$net.support,predict(loess(d.de$policy.change~d.de$net.support)), col="grey", cex=0.8)
lines(temp2[,1], temp2[,2], type="l", lty=1, lwd=3.5)
text(x=0.4, y=0.1, "UK", col='lightblue')
text(x=-0.5, y=0.6, "Germany", col='grey')
text(x=0.5, y=0.7, "Denmark", col='coral')
text(x=0.15, y=0.3, "All cases", col='black')
#dev.off()

# Logistic regression models and interaction plots (Table 1, Table A4, Figures A5:A7) --------------------------------------------
#scale and transform the variables
g$po.support.c<-g$po.support - 0.5
d$po.support.c<-d$po.support - 0.5
g$po.support.alt.c<-g$po.support.alt - 0.5
g$salience.init.log<-log(g$salience.init+1)
#d$salience.init.log<-log(d$salience.init+1)
g$country<-as.factor(g$country)
g$govpos.main.c<-g$govpos.main - 5
#g$po.polar<-ifelse(g$po.support.c>0.1 | g$po.support.c<(-0.1),1,0 )
#g$po.polar.cat<-ifelse(g$po.support.c>0.1, "2.high", ifelse(g$po.support.c<(-0.1),"0.low","1.medium"))

g.de<-subset(g, country=="DE")
#Model 1. No interactions: H1 and H2
summary(m1<-glm(event~po.support.c+gov.support+n.parties+salience.init.log+proposal+remaining+country, data=g, family='binomial'))
cbind(round(coef(m1),2),paste0("(",round(summary(m1)$coefficients[,2],2), ")"),round(summary(m1)$coefficients[,4],3))

#Model 2. Interaction between public support and government support, H3
summary(m2<-glm(event~po.support.c*gov.support+n.parties+salience.init.log+proposal+remaining+country, data=g, family='binomial'))
cbind(round(coef(m2),2),paste0("(",round(summary(m2)$coefficients[,2],2), ")"),round(summary(m2)$coefficients[,4],3))

#plot the interaction (Figure A5)
#pdf("./plots/pdf/po_gov_ie4.pdf",width=9,height=6)
#png("./plots/png/po_gov_ie4.png",width=1100, height=800, res=120)
par(mfrow=c(1,2))
e<-effect("po.support.c:gov.support", m2, se=T,confidence.level=0.75, 
          xlevels=list(po.support.c=seq(-0.5, 0.5, 0.1),gov.support=c(-2.3, 2.6) ))
plot1<-plot(e,multiline=TRUE, ylab="Probability(policy adopted)", rug=TRUE, ci="bands", xlab="Public support (centered)", main="",
            key.args=list(x=0.2, y=0.3, text=list(c("No government support", "Government support")), title="", border=F))
plot1

e2<-effect("po.support.c:gov.support", m2, se=T,confidence.level=0.75, 
           xlevels=list(po.support.c=c(-0.45, 0.45),gov.support=seq(-2.3,2.6, 0.1) ))
plot2<-plot(e2, multiline=TRUE, ylab="Probability(policy adopted)", rug=TRUE, ci="bands", xlab="Government support", main="",
            key.args=list(x=0.4, y=0.3, text=list(c("No public support", "Public support")), title="", border=F))

grid.arrange(plot1,plot2, ncol=2)
#dev.off()

#Model 3. Interaction between public support and number of parties, H5
summary(m3<-glm(event~po.support.c*n.parties+gov.support+salience.init.log+proposal+remaining+country, data=g, family='binomial'))
cbind(round(coef(m3),2),paste0("(",round(summary(m3)$coefficients[,2],2), ")"),round(summary(m3)$coefficients[,4],3))

#pdf("./plots/pdf/po_np_ie.pdf",width=9,height=6)
#png("./plots/png/po_np_ie.png",width=1100, height=800, res=120)
par(mfrow=c(1,1))
e<-effect("po.support.c:n.parties", m3, se=T,confidence.level=0.75, 
          xlevels=list(po.support.c=seq(-0.5, 0.5, 0.1),n.parties=c(1, 3) ))
plot(e,multiline=TRUE, ylab="Probability(policy adopted)", rug=TRUE, ci="bands", xlab="Public support (centered)", main="",
     key.args=list(x=0.5, y=0.2, text=list(c("Single party", "Three parties")), title="", border=F))
#dev.off()

#Model 4. Interaction between public support and salience, H4
summary(m4<-glm(event~po.support.c*salience.init.log+gov.support+n.parties+proposal+remaining+country, data=g, family='binomial'))
cbind(round(coef(m4),2),paste0("(",round(summary(m4)$coefficients[,2],2), ")"),round(summary(m4)$coefficients[,4],3))

#pdf("./plots/pdf/po_sal_ie.pdf",width=9,height=6)
#png("./plots/png/po_sal_ie.png",width=1100, height=800, res=120)
par(mfrow=c(1,1))
e<-effect("po.support.c:salience.init.log", m4, se=T,confidence.level=0.75, 
          xlevels=list(po.support.c=seq(-0.5, 0.5, 0.1),salience.init.log=c(0, 5) ))
plot(e,multiline=TRUE, ylab="Probability(policy adopted)", rug=TRUE, ci="bands", xlab="Public support (centered)", main="",
     key.args=list(x=0.5, y=0.2, text=list(c("No salience", "High salience")), title="", border=F))
#dev.off()

#Table 2, Model 5
#summary(m1de<-glm(event~po.support.c+gov.support+n.parties+salience.init.log+proposal+remaining, data=g.de, family='binomial')) #replicate the model on the german data only with Chapel Hill measures
summary(m1de<-glm(event~po.support.c+real_mean_g+n.parties+salience.init.log+proposal+remaining, data=g.de, family='binomial')) #replicate the model on the german data only with the hand-coded positions
cbind(round(coef(m1de),2),paste0("(",round(summary(m1de)$coefficients[,2],2), ")"),round(summary(m1de)$coefficients[,4],3))
#Table 2, Model 6
#summary(m2de<-glm(event~po.support.c*gov.support+n.parties+salience.init.log+proposal+remaining, data=g.de, family='binomial')) #replicate the model on the german data only with Chapel Hill measures
summary(m2de<-glm(event~po.support.c*real_mean_g+n.parties+salience.init.log+proposal+remaining, data=g.de, family='binomial')) #replicate the model on the german data only with the hand-coded positions
cbind(round(coef(m2de),2),paste0("(",round(summary(m2de)$coefficients[,2],2), ")"),round(summary(m2de)$coefficients[,4],3))

#Model A1. Model with veto player measure
summary(m5<-glm(event~po.support.c*n.parties+po.support.c*gov.vetosupport+salience.init.log+proposal+remaining+country, data=g, family='binomial'))
cbind(round(coef(m5),2),paste0("(",round(summary(m5)$coefficients[,2],2), ")"),round(summary(m5)$coefficients[,4],3))

#Model A2. Model with alternative Danish measures
summary(m6<-glm(event~po.support.c*n.parties.dk+po.support.c*gov.support.dk+salience.init.log+proposal+remaining+country, data=g, family='binomial'))
cbind(round(coef(m6),2),paste0("(",round(summary(m6)$coefficients[,2],2), ")"),round(summary(m6)$coefficients[,4],3))

#Model A3. Model with alternative PO measure
summary(m7<-glm(event~po.support.alt.c*n.parties+po.support.alt.c*gov.support+salience.init.log+proposal+remaining+country, data=g, family='binomial'))
cbind(round(coef(m7),2),paste0("(",round(summary(m7)$coefficients[,2],2), ")"),round(summary(m7)$coefficients[,4],3))

#Model A4. Full model with all three interactions
summary(m8<-glm(event~po.support.c*n.parties+po.support.c*gov.support+po.support.c*salience.init.log+proposal+remaining+country, data=g, family='binomial'))
cbind(round(coef(m8),2),paste0("(",round(summary(m8)$coefficients[,2],2), ")"),round(summary(m8)$coefficients[,4],3))


# Validation of government preferences measures (Table A2 and A3) ---------------------------
poscomp<-read.table("./data/position_comparisons.txt", sep='\t', head=T) #read the data

# Table A2.a
table(poscomp$party_position_oppose, poscomp$ch_position_final_dummy, poscomp$policy.direction) #dummies
round(prop.table(table(poscomp$party_position_oppose, poscomp$ch_position_final_dummy, poscomp$policy.direction),3),2) 

# Table A2.b
table(poscomp$party_position_oppose, poscomp$cmp_position_final_dummy, poscomp$policy.direction) #dummies
round(prop.table(table(poscomp$party_position_oppose, poscomp$cmp_position_final_dummy, poscomp$policy.direction),3),2) 

#Run regressions (Table A3)
poscomp.l<-subset(poscomp, poscomp$policy.direction==0)
poscomp.r<-subset(poscomp, poscomp$policy.direction==1)
summary(glm(party_opposition~ch_position_final, data=poscomp.l, family='binomial')) #left changes, CH data
summary(glm(party_opposition~ch_position_final, data=poscomp.r, family='binomial')) #left changes, MP data

summary(glm(party_opposition~cmp_position_final, data=poscomp.l, family='binomial')) #right changes, CH data
summary(glm(party_opposition~cmp_position_final, data=poscomp.r, family='binomial')) #right changes, MP data


# Additional varialbes for the survival analysis --------------------------
t$start<-t$month.c-1
for (i in 1:nrow(t)){
  t$gov.month[i]<-govpositions[govpositions$country==t$country[i] & govpositions$govindex==t$government[i],"edate.m"]
}

t$gov.month.c<-round((as.yearmon(t[,"month.f"], "%b %Y")-as.yearmon(t[,"gov.month"], "%b %Y"))*12, digits=0)
t$gov.sem1<-ifelse(t$gov.month.c<7,1,0)
t$gov.year1<-ifelse(t$gov.month.c<13,1,0)
t$gov.year<-factor(ifelse(t$gov.month.c<13,1,ifelse(t$gov.month.c<25,2,ifelse(t$gov.month.c<37,3, ifelse(t$gov.month.c<49,4,5)))))
t$gov.year.last<-ifelse(t$gov.month.c>35,1,0)
t$po.support.center<-t$po.support-0.5
t$po.support.alt.center<-t$po.support.alt-0.5

d$surv<-Surv(d$policy.change.duration.m.c, d$censor==0)  
d$gov.support.dummy<-ifelse(d$gov.support>0,1,0)
d$salience.init.dummy<-ifelse(d$salience.init>4,1,0)
d$govpos.main.dummy<-ifelse(d$govpos.main>5,1,0)
d$range.main.dummy<-ifelse(d$range.main>1.3,1,0)

# Survival analysis: Graphical (Figure A8) --------------------------------------------
#pdf("./plots/pdf/KMcurves.pdf",width=8,height=12)
#png("./plots/png/KMcurves.png",width=1300, height=1800, res=170)
#with public support
par(mfrow=c(3,2))
par(mar=c(5,5,2,1))

plot(survfit(surv ~ po.support.dummy, data = d, conf.type = "log-log"), lty=c(1,2), col=c("red", "blue"), xlab="Months after public opinion poll")
abline(h=seq(0,1,0.1),col="lightgrey")
abline(v=seq(0,50,10),col="lightgrey")
lines(survfit(surv ~ po.support.dummy, data = d, conf.type = "log-log"),lwd=2, lty=c(1,2), col=c("red", "blue"))
text(x=40, y=0.80,"Majority opposition", col="red")
text(x=40, y=0.5,"Majority support", col="blue")

plot(survfit(surv ~ govpos.main.dummy, data = d, conf.type = "log-log"),lty=c(1,2), col=c("red", "blue"), xlab="Months after public opinion poll")
abline(h=seq(0,1,0.1),col="lightgrey")
abline(v=seq(0,50,10),col="lightgrey")
lines(survfit(surv ~ govpos.main.dummy, data = d, conf.type = "log-log"), lty=c(1,2),lwd=2, col=c("red", "blue"))
text(x=35, y=0.80,"Right governments", col="blue")
text(x=35, y=0.48,"Left governments", col="red")

plot(survfit(surv ~ gov.support.dummy, data = d, conf.type = "log-log"),lty=c(1,2), col=c("red", "blue"), xlab="Months after public opinion poll")
abline(h=seq(0,1,0.1),col="lightgrey")
abline(v=seq(0,50,10),col="lightgrey")
lines(survfit(surv ~ gov.support.dummy, data = d, conf.type = "log-log"), lty=c(1,2),lwd=2, col=c("red", "blue"))
text(x=35, y=0.80,"Government support", col="blue")
text(x=35, y=0.48,"No government support", col="red")

plot(survfit(surv ~ range.main.dummy, data = d, conf.type = "log-log"), lty=c(1,2), col=c("red", "blue"), xlab="Months after public opinion poll")
abline(h=seq(0,1,0.1),col="lightgrey")
abline(v=seq(0,50,10),col="lightgrey")
lines(survfit(surv ~ range.main.dummy, data = d, conf.type = "log-log"), lty=c(1,2), col=c("red", "blue"), lwd=2)
text(x=41, y=0.4,"Coalition conflict low", col="red")
text(x=40, y=0.9,"Coalition conflict high", col="blue")

plot(survfit(surv ~ proposal, data = d, conf.type = "log-log"), lty=c(1,2), col=c("red", "blue"), xlab="Months after public opinion poll")
abline(h=seq(0,1,0.1),col="lightgrey")
abline(v=seq(0,50,10),col="lightgrey")
lines(survfit(surv ~ proposal, data = d, conf.type = "log-log"), lty=c(1,2), col=c("red", "blue"), lwd=2)
text(x=40, y=0.80,"No existing proposal", col="red")
text(x=40, y=0.5,"Existing proposal", col="blue")

plot(survfit(surv ~ salience.init.dummy, data = d, conf.type = "log-log"), lty=c(1,2), col=c("red", "blue"), xlab="Months after public opinion poll")
abline(h=seq(0,1,0.1),col="lightgrey")
abline(v=seq(0,50,10),col="lightgrey")
lines(survfit(surv ~ salience.init.dummy, data = d, conf.type = "log-log"), lty=c(1,2), col=c("red", "blue"), lwd=2)
text(x=40, y=0.80,"Low salience", col="red")
text(x=40, y=0.5,"High salience", col="blue")
#dev.off()
# Survival analysis models with time-varying covaraites (Table A5)----------------------------------------------
summary(ma5<-coxph(Surv(start, month.c, event)~po.support.center+coal.conflict.main+gov.support.t+
                    log(salience.init+1)+proposal+
                    strata(country)+cluster(case.n), data=t))
round(summary(ma5)$coef[,c(1,4,6)], digits=2)


summary(ma6<-coxph(Surv(start, month.c, event)~po.support.center*coal.conflict.main+gov.support.t+
                    log(salience.init+1)+proposal+
                    strata(country)+cluster(case.n), data=t))
round(summary(ma6)$coef[,c(1,4,6)], digits=2)

summary(ma7<-coxph(Surv(start, month.c, event)~po.support.center+coal.conflict.main+po.support.center*gov.support.t+
                    log(salience.init+1)+proposal+
                    strata(country)+cluster(case.n), data=t))
round(summary(ma7)$coef[,c(1,4,6)], digits=2)
###THE END